#BiocManager::install("clusterProfiler")
#install.packages("enrichplot")
#installed.packages("patchwork")
#BiocManager::install("SingleR")
#BiocManager::install("celldex")
#BiocManager::install("ensembldb")
#BiocManager::install('limma')
#BiocManager::install("DESeq2")
#organism = "org.Mm.eg.db"
#BiocManager::install(organism)
#BiocManager::install("pathview")
#install.packages("Seurat")
#install.packages("remotes")
#remotes::install_github("satijalab/seurat")
#BiocManager::install("rhdf5")
#install.packages("hdf5r")
#BiocManager::install("NeuCA")
#install.packages('devtools')
#install_github("ggjlab/scMCA")
#install.packages("pheatmap")
#BiocManager::install("scater")
#devtools::install_github("sqjin/CellChat")
#BiocManager::install("ComplexHeatmap")
#BiocManager::install("SingleCellExperiment")
In this step, we are loading the libraries which are required for analysis and plotting of the data.
library(Seurat)
library(devtools)
library(ggplot2)
library(patchwork)
library(dplyr)
library(clusterProfiler)
library(enrichplot)
# we use ggplot2 to add x axis labels (ex: ridgeplot)
library(stringr)
library(celldex)
library(ensembldb)
library(SingleR)
library(DESeq2)
library(org.Mm.eg.db)
library(rhdf5)
library(hdf5r)
library(pheatmap)
library(ComplexHeatmap)
library(CellChat)
library(patchwork)
options(stringsAsFactors = FALSE)
library(NMF)
library(ggalluvial)
my_cols <- c('B cells'='#F8766D','DC'='#39568CFF','Endothelial cells'='#CD9600','Epithelial cells'='#00C19A','Fibroblasts'= "#C77CFF",
'Macrophages'='#00A9FF','Monocytes'='#0CB702','Stromal cells'='#FF61CC', "ILC" = "grey", "Neutrophils" = "Darkred")
# Mouse 1
ntm_mouse_1 <- "/Volumes/Seagate_4TB/spatial_ntm_samples/demultiplexed/ntm_sample_round_one/outs/"
list.files(ntm_mouse_1)
## [1] "analysis" "cloupe.cloupe"
## [3] "filtered_feature_bc_matrix_ntm_1" "filtered_feature_bc_matrix.h5"
## [5] "metrics_summary.csv" "molecule_info.h5"
## [7] "possorted_genome_bam.bam" "possorted_genome_bam.bam.bai"
## [9] "probe_set.csv" "raw_feature_bc_matrix"
## [11] "raw_feature_bc_matrix.h5" "spatial"
## [13] "spatial_enrichment.csv" "web_summary.html"
ntm_mouse_1_data <- Load10X_Spatial(
ntm_mouse_1,
filename = "filtered_feature_bc_matrix.h5",
assay = "spatial",
slice = "ntm",
filter.matrix = TRUE,
to.upper = FALSE,
image = NULL)
# Mouse 2
ntm_mouse_2 <- "/Volumes/Seagate_4TB/spatial_ntm_samples/demultiplexed/ntm_sample_round_two/outs/"
list.files(ntm_mouse_2)
## [1] "analysis" "cloupe.cloupe"
## [3] "filtered_feature_bc_matrix_ntm_2" "filtered_feature_bc_matrix.h5"
## [5] "metrics_summary.csv" "molecule_info.h5"
## [7] "ntm_only_web_summary.html" "possorted_genome_bam.bam"
## [9] "possorted_genome_bam.bam.bai" "probe_set.csv"
## [11] "raw_feature_bc_matrix" "raw_feature_bc_matrix.h5"
## [13] "spatial" "spatial_enrichment.csv"
ntm_mouse_2_data <- Load10X_Spatial(
ntm_mouse_2,
filename = "filtered_feature_bc_matrix.h5",
assay = "spatial",
slice = "ntm",
filter.matrix = TRUE,
to.upper = FALSE,
image = NULL)
# Mouse 1
ntm_mouse_1_data <- SCTransform(ntm_mouse_1_data, assay = "spatial", verbose = FALSE)
# Mouse 2
ntm_mouse_2_data <- SCTransform(ntm_mouse_2_data, assay = "spatial", verbose = FALSE)
p1 <- SpatialFeaturePlot(ntm_mouse_1_data, features = "Cd19", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p2 <- SpatialFeaturePlot(ntm_mouse_1_data, features = "Igha", alpha = c(0.1, 1), pt.size.factor = 1.5)
p3 <- SpatialFeaturePlot(ntm_mouse_1_data, features = "Bcl6", alpha = c(0.1, 1), pt.size.factor = 1.5)
p4 <- SpatialFeaturePlot(ntm_mouse_1_data, features = "Lta", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p5 <- SpatialFeaturePlot(ntm_mouse_1_data, features = "Ltb", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p6 <- SpatialFeaturePlot(ntm_mouse_1_data, features = "Cd4", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p7 <- SpatialFeaturePlot(ntm_mouse_1_data, features = "Mki67", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p8 <- SpatialFeaturePlot(ntm_mouse_1_data, features = "Cxcr5", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8
p1 <- SpatialFeaturePlot(ntm_mouse_2_data, features = "Cd19", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p2 <- SpatialFeaturePlot(ntm_mouse_2_data, features = "Igha", alpha = c(0.1, 1), pt.size.factor = 1.5)
p3 <- SpatialFeaturePlot(ntm_mouse_2_data, features = "Bcl6", alpha = c(0.1, 1), pt.size.factor = 1.5)
p4 <- SpatialFeaturePlot(ntm_mouse_2_data, features = "Lta", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p5 <- SpatialFeaturePlot(ntm_mouse_2_data, features = "Ltb", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p6 <- SpatialFeaturePlot(ntm_mouse_2_data, features = "Cd4", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p7 <- SpatialFeaturePlot(ntm_mouse_2_data, features = "Mki67", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p8 <- SpatialFeaturePlot(ntm_mouse_2_data, features = "Cxcr5", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8
# Mouse 1
ntm_mouse_1_analyzed <- RunPCA(ntm_mouse_1_data, assay = "SCT", verbose = FALSE)
ntm_mouse_1_analyzed <- FindNeighbors(ntm_mouse_1_analyzed, reduction = "pca", dims = 1:5)
ntm_mouse_1_analyzed <- FindClusters(ntm_mouse_1_analyzed, verbose = FALSE, resolution = 0.8)
ntm_mouse_1_analyzed <- RunUMAP(ntm_mouse_1_analyzed, reduction = "pca", dims = 1:5)
# Mouse 2
ntm_mouse_2_analyzed <- RunPCA(ntm_mouse_2_data, assay = "SCT", verbose = FALSE)
ntm_mouse_2_analyzed <- FindNeighbors(ntm_mouse_2_analyzed, reduction = "pca", dims = 1:5)
ntm_mouse_2_analyzed <- FindClusters(ntm_mouse_2_analyzed, verbose = FALSE, resolution = 0.8)
ntm_mouse_2_analyzed <- RunUMAP(ntm_mouse_2_analyzed, reduction = "pca", dims = 1:5)
# First convert the seurat object to single cell experiment object (otherwise the SingleR pipeline will not work)
# Mouse 1
ntm_mouse_1_analyzed.sce <- as.SingleCellExperiment(ntm_mouse_1_analyzed)
# Mouse 2
ntm_mouse_2_analyzed.sce <- as.SingleCellExperiment(ntm_mouse_2_analyzed)
#create reference data
ref.data <- celldex::ImmGenData()
# run singleR pipeline to find the cell types
## Mouse 1
cell_types_mouse_1 <- SingleR(test=ntm_mouse_1_analyzed.sce, assay.type.test=1,
ref=ref.data, labels=ref.data$label.main)
## Mouse 2
cell_types_mouse_2 <- SingleR(test=ntm_mouse_2_analyzed.sce, assay.type.test=1,
ref=ref.data, labels=ref.data$label.main)
# view cell types
# Mouse 1
table(cell_types_mouse_1$labels)
##
## B cells DC Endothelial cells Epithelial cells
## 152 59 868 35
## Fibroblasts Macrophages Monocytes Neutrophils
## 716 476 7 1
## Stromal cells
## 601
#checking cell types Mouse 1
plotScoreHeatmap(cell_types_mouse_1)
plotDeltaDistribution(cell_types_mouse_1)
summary(is.na(cell_types_mouse_1$pruned.labels))
## Mode FALSE TRUE
## logical 2901 14
# Mouse 2
table(cell_types_mouse_2$labels)
##
## B cells DC Endothelial cells Epithelial cells
## 124 19 970 28
## Fibroblasts Macrophages Monocytes Neutrophils
## 234 477 6 4
## Stromal cells
## 312
#checking cell types Mouse 2
plotScoreHeatmap(cell_types_mouse_2)
plotDeltaDistribution(cell_types_mouse_2)
summary(is.na(cell_types_mouse_2$pruned.labels))
## Mode FALSE TRUE
## logical 2167 7
# Mouse 1
ntm_mouse_1_analyzed[["ref.data"]] <- cell_types_mouse_1$labels
# Mouse 2
ntm_mouse_2_analyzed[["ref.data"]] <- cell_types_mouse_2$labels
p1 <- DimPlot(ntm_mouse_1_analyzed, group.by = c("ref.data"), reduction = "umap",label = F, pt.size = 1.5, label.size = 5, cols = my_cols)
p2 <- SpatialDimPlot(ntm_mouse_1_analyzed, group.by = c("ref.data"), label = FALSE, label.size = 5, cols = my_cols)
p1 + p2
p1 <- DimPlot(ntm_mouse_2_analyzed, group.by = c("ref.data"), reduction = "umap",label = F, pt.size = 1.5, label.size = 5, cols = my_cols)
p2 <- SpatialDimPlot(ntm_mouse_2_analyzed, group.by = c("ref.data"), label = FALSE, label.size = 5, cols = my_cols)
p1 + p2
VlnPlot(ntm_mouse_1_analyzed, group.by = c("ref.data"), features = c("Cxcr5", "Mki67", "Ltb", "Cd38", "Ccr6", "Tnfrsf13c"), cols = my_cols)
VlnPlot(ntm_mouse_2_analyzed, group.by = c("ref.data"), features = c("Cxcr5", "Mki67", "Ltb", "Cd38", "Ccr6", "Tnfrsf13c"), cols = my_cols)
B_cell_mouse_1 <- FindMarkers(ntm_mouse_1_analyzed, group.by = "ref.data", ident.1 = "B cells", min.pct = 0.25)
head(B_cell_mouse_1, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Tnfrsf13c 4.554733e-117 1.255106 0.816 0.142 8.323320e-113
## Ms4a1 3.671873e-111 1.080410 0.724 0.107 6.709980e-107
## Cd79b 2.869009e-105 1.819936 0.914 0.254 5.242828e-101
## Iglc3 9.359354e-105 1.368350 0.868 0.201 1.710328e-100
## Pou2af1 1.805543e-101 1.920921 0.974 0.349 3.299449e-97
## Ptprcap 1.197781e-99 1.093994 0.822 0.170 2.188825e-95
## Cd79a 2.339769e-98 2.275292 0.987 0.425 4.275694e-94
## Ltb 7.343352e-95 1.289891 0.875 0.220 1.341924e-90
## Spib 9.978568e-94 1.299672 0.763 0.157 1.823483e-89
## Cd19 2.272512e-93 1.590420 0.849 0.222 4.152788e-89
write.csv(B_cell_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/B_cell_mouse_1.csv")
DC_mouse_1 <- FindMarkers(ntm_mouse_1_analyzed, group.by = "ref.data", ident.1 = "DC", min.pct = 0.25)
head(DC_mouse_1, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## H2-Eb1 5.334939e-24 1.0919341 1.000 0.994 9.749067e-20
## C3 4.157713e-23 1.0263028 1.000 0.961 7.597804e-19
## Tmsb4x 1.105303e-21 0.7047462 1.000 1.000 2.019830e-17
## H2-Aa 1.637641e-21 1.0455447 1.000 0.887 2.992626e-17
## Psmb8 4.052254e-21 0.8606870 1.000 0.818 7.405089e-17
## Ccl22 1.634381e-20 0.6154602 0.525 0.123 2.986668e-16
## Inmt 2.315338e-20 -1.8594434 0.458 0.870 4.231048e-16
## Cd52 7.290028e-20 0.9789633 0.983 0.765 1.332180e-15
## H2-Ab1 8.191305e-20 0.8728392 1.000 0.960 1.496879e-15
## Serpina3g 1.069894e-19 1.0868829 1.000 0.829 1.955125e-15
write.csv(DC_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/DC_mouse_1.csv")
Stromal_cells_mouse_1 <- FindMarkers(ntm_mouse_1_analyzed, group.by = "ref.data", ident.1 = "Stromal cells", min.pct = 0.25)
head(Stromal_cells_mouse_1, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Myh11 6.158379e-36 0.5558271 0.532 0.285 1.125382e-31
## Gsn 3.134160e-33 0.3994196 0.993 0.925 5.727364e-29
## Des 1.691635e-32 0.5752681 0.654 0.432 3.091294e-28
## Ctss 2.689407e-32 -0.6236161 0.952 0.983 4.914621e-28
## Tagln 2.766731e-28 0.4075501 0.461 0.249 5.055924e-24
## Myl9 5.016710e-27 0.4943200 0.496 0.284 9.167536e-23
## Car3 1.236110e-25 1.1447737 0.374 0.189 2.258867e-21
## B2m 3.247562e-25 -0.3333778 0.998 1.000 5.934595e-21
## Tmsb4x 1.472912e-24 -0.3007193 1.000 1.000 2.691599e-20
## Synpo2 1.333461e-23 0.3428623 0.381 0.197 2.436767e-19
write.csv(Stromal_cells_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Stromal_cells_mouse_1.csv")
Endothelial_cells_mouse_1 <- FindMarkers(ntm_mouse_1_analyzed, group.by = "ref.data", ident.1 = "Endothelial cells", min.pct = 0.25)
head(Endothelial_cells_mouse_1, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Cldn5 1.823586e-102 0.6927932 0.992 0.812 3.332422e-98
## Igfbp2 3.821331e-95 0.7524716 0.983 0.778 6.983101e-91
## Epas1 8.978855e-91 0.6125308 0.995 0.899 1.640796e-86
## Ager 1.241455e-86 0.5622126 1.000 0.936 2.268635e-82
## Acer2 8.633030e-86 0.6179743 0.983 0.809 1.577600e-81
## Clic5 2.306495e-85 0.5693432 0.993 0.888 4.214889e-81
## Egfl7 2.202114e-82 0.5948463 0.984 0.834 4.024143e-78
## C3 9.870730e-82 -0.8186092 0.937 0.972 1.803777e-77
## Vegfa 8.421764e-80 0.5778378 0.988 0.865 1.538993e-75
## Ace 3.230574e-78 0.6225789 0.963 0.766 5.903551e-74
write.csv(Endothelial_cells_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Endothelial_cells_mouse_1.csv")
Epithelial_cell_mouse_1 <- FindMarkers(ntm_mouse_1_analyzed, group.by = "ref.data", ident.1 = "Epithelial cells", min.pct = 0.25)
head(Epithelial_cell_mouse_1, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Gp2 3.057056e-67 1.0931770 0.686 0.044 5.586464e-63
## Tmc5 4.457791e-65 0.7977855 0.600 0.034 8.146167e-61
## Cd177 2.446346e-56 0.8906558 0.629 0.044 4.470452e-52
## B3galt2 1.335568e-51 0.5792686 0.514 0.031 2.440617e-47
## Bpifb1 3.657771e-51 1.9902073 0.800 0.088 6.684210e-47
## Erich3 3.941690e-50 0.9298852 0.800 0.082 7.203045e-46
## Muc5b 1.797909e-47 1.3634905 0.686 0.066 3.285499e-43
## Pklr 4.992440e-47 0.8550121 0.629 0.053 9.123184e-43
## Tff2 1.128400e-43 1.4010201 0.686 0.072 2.062038e-39
## Kcnk12 1.811543e-42 0.7501989 0.657 0.064 3.310413e-38
write.csv(Epithelial_cell_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Epithelial_cell_mouse_1.csv")
Fibroblasts_mouse_1 <- FindMarkers(ntm_mouse_1_analyzed, group.by = "ref.data", ident.1 = "Fibroblasts", min.pct = 0.25)
head(Fibroblasts_mouse_1, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Gsn 2.679271e-39 0.4309985 0.997 0.920 4.896101e-35
## Psap 8.682005e-31 -0.3999956 0.999 1.000 1.586550e-26
## Inmt 1.692121e-30 0.4104036 0.980 0.824 3.092182e-26
## Ccl5 9.335028e-27 -0.6195179 0.538 0.692 1.705883e-22
## Ltbp4 8.209228e-25 0.3411285 0.939 0.819 1.500154e-20
## Ogn 1.929405e-24 0.3335998 0.747 0.543 3.525795e-20
## Ctss 1.099770e-23 -0.5534927 0.973 0.978 2.009719e-19
## Socs1 3.691625e-23 -0.4632449 0.409 0.584 6.746076e-19
## H2-Ab1 2.379039e-22 -0.4181848 0.946 0.965 4.347456e-18
## Fmo2 2.668078e-22 0.3470290 0.957 0.808 4.875645e-18
write.csv(Fibroblasts_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Fibroblasts_mouse_1.csv")
Macrophages_mouse_1 <- FindMarkers(ntm_mouse_1_analyzed, group.by = "ref.data", ident.1 = "Macrophages", min.pct = 0.25)
head(Macrophages_mouse_1, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Ctss 4.554123e-140 1.2138975 1.000 0.972 8.322204e-136
## Inmt 9.007957e-120 -1.5120459 0.620 0.909 1.646114e-115
## Psap 2.867758e-116 0.8116796 1.000 0.999 5.240540e-112
## B2m 2.043846e-113 0.6627254 1.000 1.000 3.734924e-109
## Gsn 4.765322e-111 -1.0967002 0.836 0.959 8.708150e-107
## Gbp2b 4.006186e-110 1.1110715 0.979 0.754 7.320904e-106
## Tnfaip2 5.875226e-108 1.1302766 0.968 0.771 1.073639e-103
## H2-Eb1 3.943305e-107 0.7819851 1.000 0.993 7.205996e-103
## H2-Ab1 1.595128e-106 0.8231051 1.000 0.953 2.914937e-102
## Cd74 1.150780e-104 0.5097265 1.000 1.000 2.102936e-100
write.csv(Macrophages_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Macrophages_mouse_1.csv")
Monocytes_mouse_1 <- FindMarkers(ntm_mouse_1_analyzed, group.by = "ref.data", ident.1 = "Monocytes", min.pct = 0.25)
head(Monocytes_mouse_1, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Nlgn3 1.740289e-15 0.3512044 0.286 0.008 3.180203e-11
## Olfr447 6.404532e-15 0.3507122 0.286 0.008 1.170364e-10
## Jph3 4.958409e-13 0.3487454 0.286 0.010 9.060997e-09
## Cma1 6.626389e-12 0.3467812 0.286 0.011 1.210906e-07
## Polr3f 1.413811e-11 0.3462906 0.286 0.011 2.583598e-07
## Zfp846 2.843161e-11 0.3462906 0.286 0.011 5.195592e-07
## Zfp982 3.118175e-11 0.4783202 0.429 0.025 5.698154e-07
## Gm15056 3.191898e-11 0.8252492 0.571 0.047 5.832874e-07
## U2af1l4 8.187911e-11 0.4754203 0.429 0.026 1.496259e-06
## Adgrg3 1.059973e-10 0.3453098 0.286 0.012 1.936995e-06
write.csv(Monocytes_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Monocytes_mouse_1.csv")
B_cell_mouse_2 <- FindMarkers(ntm_mouse_2_analyzed, group.by = "ref.data", ident.1 = "B cells", min.pct = 0.25)
head(B_cell_mouse_2, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Bank1 4.202720e-83 1.3287032 0.726 0.127 5.797232e-79
## Coch 1.488250e-82 0.4608521 0.274 0.008 2.052892e-78
## Tnfrsf13c 1.201551e-81 1.7538640 0.855 0.208 1.657419e-77
## Cxcr5 2.614613e-81 1.0992235 0.653 0.096 3.606597e-77
## Fcrl5 1.683029e-79 0.8158863 0.508 0.051 2.321570e-75
## Ccr6 4.425674e-78 1.5810476 0.790 0.177 6.104775e-74
## Cd79a 1.041844e-77 2.5173171 1.000 0.572 1.437120e-73
## Cd79b 1.195236e-77 2.1139509 0.984 0.383 1.648708e-73
## Cd19 1.761850e-77 1.8569971 0.831 0.207 2.430296e-73
## Spib 2.477439e-76 1.7252838 0.855 0.229 3.417379e-72
write.csv(B_cell_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/B_cell_mouse_2.csv")
DC_mouse_2 <- FindMarkers(ntm_mouse_2_analyzed, group.by = "ref.data", ident.1 = "DC", min.pct = 0.25)
head(DC_mouse_2, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Xcr1 3.062816e-11 0.3883528 0.368 0.044 4.224848e-07
## Foxd2 8.827878e-11 0.3472199 0.316 0.034 1.217717e-06
## Ccr7 1.753044e-10 0.7015559 0.789 0.206 2.418149e-06
## Ints10 2.099668e-10 0.3006765 0.263 0.025 2.896283e-06
## Ttc39a 3.705338e-10 0.5240484 0.684 0.149 5.111144e-06
## Mettl25 4.552337e-10 0.4161537 0.421 0.063 6.279494e-06
## Cabp4 1.689226e-09 0.3401178 0.316 0.039 2.330119e-05
## Nckap1l 3.071381e-09 0.4610973 0.421 0.070 4.236663e-05
## Adam11 6.130534e-09 0.3530569 0.263 0.030 8.456459e-05
## Cd5 7.717586e-09 0.2941632 0.263 0.030 1.064564e-04
write.csv(DC_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/DC_mouse_2.csv")
Stromal_cells_mouse_2 <- FindMarkers(ntm_mouse_2_analyzed, group.by = "ref.data", ident.1 = "Stromal cells", min.pct = 0.25)
head(Stromal_cells_mouse_2, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Myh11 8.950958e-39 0.9616085 0.657 0.331 1.234695e-34
## Tagln 5.600725e-32 0.9005221 0.619 0.328 7.725641e-28
## Acta2 6.307122e-32 1.0824300 0.679 0.385 8.700044e-28
## Ltbp4 6.419828e-30 0.5857921 0.952 0.800 8.855510e-26
## Myl9 1.246700e-27 0.7490902 0.577 0.302 1.719698e-23
## Ptgis 2.324078e-26 0.6232297 0.881 0.669 3.205833e-22
## Eln 8.537455e-25 0.7410893 0.904 0.789 1.177657e-20
## Ppp1r14a 4.089734e-23 0.4813196 0.990 0.826 5.641379e-19
## Gsn 3.191197e-22 0.4398570 0.994 0.899 4.401937e-18
## Ccn2 3.356066e-22 0.6573885 0.853 0.664 4.629358e-18
write.csv(Stromal_cells_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Stromal_cells_mouse_2.csv")
Endothelial_cells_mouse_2 <- FindMarkers(ntm_mouse_2_analyzed, group.by = "ref.data", ident.1 = "Endothelial cells", min.pct = 0.25)
head(Endothelial_cells_mouse_2, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Igfbp2 5.232736e-141 1.1102123 0.992 0.727 7.218035e-137
## Cldn5 3.910543e-132 0.9235680 0.996 0.804 5.394203e-128
## Epas1 2.721131e-120 0.8123111 0.999 0.890 3.753528e-116
## Ager 6.818852e-113 0.6825610 1.000 0.950 9.405924e-109
## Pakap.1 1.259013e-109 0.7385255 1.000 0.883 1.736682e-105
## Clic5 1.048858e-105 0.6633032 0.999 0.914 1.446795e-101
## Ace 4.903890e-103 0.7553659 0.982 0.741 6.764425e-99
## C3 1.469702e-98 -0.7505188 0.999 0.997 2.027308e-94
## Tns1 1.850102e-97 0.6977257 0.998 0.880 2.552031e-93
## Egfl7 1.363452e-95 0.6908519 0.995 0.827 1.880745e-91
write.csv(Endothelial_cells_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Endothelial_cells_mouse_2.csv")
Epithelial_cell_mouse_2 <- FindMarkers(ntm_mouse_2_analyzed, group.by = "ref.data", ident.1 = "Epithelial cells", min.pct = 0.25)
head(Epithelial_cell_mouse_2, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Erich3 5.496562e-20 0.4546580 0.393 0.039 7.581958e-16
## Cxcl17 2.337558e-17 0.8259062 0.571 0.100 3.224428e-13
## Tex26 2.609954e-17 0.5173289 0.393 0.045 3.600170e-13
## Bpifb1 4.200883e-17 0.5399574 0.286 0.025 5.794698e-13
## Ubxn10 4.527741e-17 0.8378696 0.643 0.127 6.245566e-13
## Lrrc71 9.566381e-17 0.5689175 0.500 0.075 1.319587e-12
## Ankfn1 3.546138e-16 0.3531888 0.321 0.032 4.891543e-12
## Zfp750 7.178780e-16 0.4780801 0.393 0.049 9.902410e-12
## Dynlrb2 1.000800e-15 0.5289973 0.429 0.059 1.380503e-11
## Tmem212 1.184842e-15 0.6631443 0.571 0.103 1.634371e-11
write.csv(Epithelial_cell_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Epithelial_cell_mouse_2.csv")
Fibroblasts_mouse_2 <- FindMarkers(ntm_mouse_2_analyzed, group.by = "ref.data", ident.1 = "Fibroblasts", min.pct = 0.25)
head(Fibroblasts_mouse_2, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Cst3 3.397128e-24 0.3258673 1.000 1.000 4.685998e-20
## Clec3b 2.544528e-20 0.5677545 0.932 0.746 3.509922e-16
## Gsn 3.803288e-19 0.4445807 0.996 0.903 5.246255e-15
## Psap 3.958205e-17 -0.4074829 1.000 1.000 5.459949e-13
## Ctss 4.568216e-14 -0.4660153 1.000 1.000 6.301398e-10
## Tppp3 5.760751e-14 0.5288847 0.991 0.909 7.946379e-10
## Cybb 1.130674e-13 -0.4894617 0.872 0.926 1.559652e-09
## Lgals3 1.663983e-13 -0.5059102 0.979 0.974 2.295298e-09
## Gbp2b 2.252567e-12 -0.5267407 0.932 0.974 3.107191e-08
## Ctsz 2.496011e-12 -0.4260957 0.927 0.968 3.442998e-08
write.csv(Fibroblasts_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Fibroblasts_mouse_2.csv")
Macrophages_mouse_2 <- FindMarkers(ntm_mouse_2_analyzed, group.by = "ref.data", ident.1 = "Macrophages", min.pct = 0.25)
head(Macrophages_mouse_2, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Igfbp2 3.409736e-105 -1.5826761 0.627 0.906 4.703390e-101
## Tnfaip2 3.424259e-100 1.1606650 1.000 0.978 4.723423e-96
## Epas1 2.429809e-99 -1.0990461 0.878 0.955 3.351679e-95
## Inmt 3.250042e-99 -1.3557853 0.610 0.906 4.483108e-95
## Gbp2b 3.785450e-99 0.9670132 0.996 0.962 5.221650e-95
## Psap 1.318899e-97 0.6908248 1.000 1.000 1.819289e-93
## Ctss 1.242439e-96 0.8073770 1.000 1.000 1.713820e-92
## Tns1 1.760472e-94 -1.0432070 0.870 0.951 2.428395e-90
## Fmo2 9.261547e-92 -1.1333143 0.618 0.902 1.277538e-87
## Cldn5 8.194109e-91 -1.1740409 0.776 0.922 1.130295e-86
write.csv(Macrophages_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Macrophages_mouse_2.csv")
Monocytes_mouse_2 <- FindMarkers(ntm_mouse_2_analyzed, group.by = "ref.data", ident.1 = "Monocytes", min.pct = 0.25)
head(Monocytes_mouse_2, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Rnf207 7.675060e-54 0.4123782 0.333 0.002 1.058698e-49
## Nup210l 7.992847e-30 0.4090608 0.333 0.004 1.102533e-25
## P2rx2 4.105985e-22 0.4064125 0.333 0.006 5.663795e-18
## Tdrkh 1.039971e-14 0.4011303 0.333 0.010 1.434535e-10
## Ssc5d 3.216632e-12 0.3978387 0.333 0.012 4.437023e-08
## Cdc6 8.956139e-12 0.5442835 0.500 0.028 1.235410e-07
## 2410131K14Rik 4.319930e-10 0.6535795 0.667 0.059 5.958912e-06
## Clstn2 3.541658e-09 0.6435651 0.667 0.065 4.885364e-05
## Crip3 9.169083e-09 0.3899695 0.333 0.018 1.264783e-04
## Isoc2b 5.081043e-08 0.3867032 0.333 0.019 7.008790e-04
write.csv(Monocytes_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Monocytes_mouse_2.csv")
Neutrophils_mouse_2 <- FindMarkers(ntm_mouse_2_analyzed, group.by = "ref.data", ident.1 = "Neutrophils", min.pct = 0.25)
head(Neutrophils_mouse_2, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Sema5b 4.443551e-25 0.3192712 0.25 0.002 6.129434e-21
## Cyp26a1 4.443551e-25 0.3192712 0.25 0.002 6.129434e-21
## Frmd7 4.443551e-25 0.3192712 0.25 0.002 6.129434e-21
## Zfp469 4.195789e-21 0.3186077 0.25 0.002 5.787671e-17
## Spock1 4.195789e-21 0.3186077 0.25 0.002 5.787671e-17
## Med12 2.364461e-18 0.5809790 0.25 0.003 3.261537e-14
## Cdx2 2.932706e-18 0.3179446 0.25 0.003 4.045375e-14
## Hoxc5 2.932706e-18 0.3179446 0.25 0.003 4.045375e-14
## Cartpt 3.039891e-18 0.3152950 0.25 0.003 4.193226e-14
## Grm5 5.413888e-18 0.5677795 0.50 0.012 7.467918e-14
write.csv(Neutrophils_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Neutrophils_mouse_2.csv")
B_cell_deseq <- FindMarkers(ntm_mouse_1_analyzed, group.by = "ref.data", ident.1 = "Macrophages", test.use = "DESeq2")
write.csv(B_cell_deseq, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/B_cell_deseq.csv")
ntm_cellchat <- createCellChat(object = ntm_mouse_2_analyzed, group.by = "ref.data")
## The `data` slot in the default assay is used. The default assay is SCT
## The `meta.data` slot in the Seurat object is used as cell meta information
## The cell groups used for CellChat analysis are B cells DC Endothelial cells Epithelial cells Fibroblasts Macrophages Monocytes Neutrophils Stromal cells
CellChatDB <- CellChatDB.mouse
showDatabaseCategory(CellChatDB)
CellChatDB.use <- CellChatDB
ntm_cellchat@DB <- CellChatDB.use
ntm_cellchat <- subsetData(ntm_cellchat)
## Issue identified!! Please check the official Gene Symbol of the following genes:
## H2-BI H2-Ea-ps
ntm_cellchat <- identifyOverExpressedGenes(ntm_cellchat)
ntm_cellchat <- identifyOverExpressedInteractions(ntm_cellchat)
ntm_cellchat <- computeCommunProb(ntm_cellchat)
ntm_cellchat <- filterCommunication(ntm_cellchat, min.cells = 10)
## The cell-cell communication related with the following cell groups are excluded due to the few number of cells: Monocytes Neutrophils
ntm_cellchat <- computeCommunProbPathway(ntm_cellchat)
ntm_cellchat <- aggregateNet(ntm_cellchat)
groupSize <- as.numeric(table(ntm_cellchat@idents))
par(mfrow = c(1,2), xpd=TRUE)
netVisual_circle(ntm_cellchat@net$count, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Number of interactions")
netVisual_circle(ntm_cellchat@net$weight, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Interaction weights/strength")
# See interaction of each cell type with the another
mat1 <- ntm_cellchat@net$weight
par(mfrow = c(3,4), xpd=TRUE)
for (i in 1:nrow(mat1)) {
mat2 <- matrix(0, nrow = nrow(mat1), ncol = ncol(mat1), dimnames = dimnames(mat1))
mat2[i, ] <- mat1[i, ]
netVisual_circle(mat2, vertex.weight = groupSize, weight.scale = T, edge.weight.max = max(mat1), title.name = rownames(mat1)[i])
}
##Visualize signalling pathways in each cell
# Identify signaling pathways showing significant communications
ntm_cellchat@netP$pathways
## [1] "COLLAGEN" "APP" "THBS" "FN1" "MIF"
## [6] "COMPLEMENT" "GALECTIN" "MHC-II" "VEGF" "LAMININ"
## [11] "SEMA4" "ICAM" "CXCL" "GAS" "CCL"
## [16] "ITGAL-ITGB2" "JAM" "PECAM1" "CD22" "CD45"
## [21] "PDGF" "CD52" "TENASCIN" "THY1" "TGFb"
## [26] "FGF" "SEMA3" "CHEMERIN" "PD-L1" "LT"
## [31] "SEMA7" "SPP1" "NOTCH" "GRN" "SEMA6"
## [36] "ESAM" "IL16" "EPHB" "PARs" "CDH"
## [41] "ncWNT" "IGF" "CD39" "APRIL" "ICOS"
## [46] "PDL2" "NRG" "XCR" "IL1"
# Evaluate one of the pathways
pathways.show <- c("MHC-II")
vertex.receiver = seq(2,5) # a numeric vector.
?netVisual_aggregate(ntm_cellchat, signaling = pathways.show, vertex.receiver = vertex.receiver)
# Chord diagram
par(mfrow=c(1,1))
netVisual_aggregate(ntm_cellchat, signaling = pathways.show, layout = "chord")
# Heatmap
par(mfrow=c(1,1))
netVisual_heatmap(ntm_cellchat, signaling = pathways.show, color.heatmap = "Reds")
netAnalysis_contribution(ntm_cellchat, signaling = pathways.show)
pairLR.MHC <- extractEnrichedLR(ntm_cellchat, signaling = pathways.show, geneLR.return = FALSE)
LR.show <- pairLR.MHC[1,] # show one ligand-receptor pair
# Hierarchy plot
vertex.receiver = seq(1,4) # a numeric vector
netVisual_individual(ntm_cellchat, signaling = pathways.show, pairLR.use = LR.show, vertex.receiver = vertex.receiver)
## [[1]]
# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
netVisual_bubble(ntm_cellchat, sources.use = 1, targets.use = c(2:8), remove.isolate = FALSE)
# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
# show all the interactions sending from Inflam.FIB
netVisual_chord_gene(ntm_cellchat, sources.use = 1, targets.use = c(2:8), lab.cex = 0.5,legend.pos.y = 30)
ntm_cellchat <- netAnalysis_computeCentrality(ntm_cellchat, slot.name = "netP")
# the slot 'netP' means the inferred intercellular communication network of signaling pathways
# Visualize the computed centrality scores using heatmap, allowing ready identification of major signaling roles of cell groups
netAnalysis_signalingRole_network(ntm_cellchat, signaling = pathways.show, width = 8, height = 2.5, font.size = 10)
ht1 <- netAnalysis_signalingRole_heatmap(ntm_cellchat, pattern = "outgoing")
ht2 <- netAnalysis_signalingRole_heatmap(ntm_cellchat, pattern = "incoming")
ht1 + ht2
selectK(ntm_cellchat, pattern = "outgoing")
nPatterns = 3
ntm_cellchat <- identifyCommunicationPatterns(ntm_cellchat, pattern = "outgoing", k = nPatterns)
# river plot
netAnalysis_river(ntm_cellchat, pattern = "outgoing")
# dot plot
netAnalysis_dot(ntm_cellchat, pattern = "outgoing")
selectK(ntm_cellchat, pattern = "incoming")
nPatterns = 3
ntm_cellchat <- identifyCommunicationPatterns(ntm_cellchat, pattern = "incoming", k = nPatterns)
# river plot
netAnalysis_river(ntm_cellchat, pattern = "incoming")
# dot plot
netAnalysis_dot(ntm_cellchat, pattern = "incoming")
organism = org.Mm.eg.db
# reading in data from deseq2
df = read.csv("/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/B_cell_deseq.csv", header=TRUE)
# we want the log2 fold change
original_gene_list <- df$avg_log2FC
# name the vector
names(original_gene_list) <- df$X
# omit any NA values
gene_list<-na.omit(original_gene_list)
# sort the list in decreasing order (required for clusterProfiler)
gene_list = sort(gene_list, decreasing = TRUE)
#gene set enrichment function for B cell cluster
gse <- gseGO(geneList=gene_list,
ont ="ALL",
keyType = "SYMBOL",
minGSSize = 3,
maxGSSize = 800,
pvalueCutoff = 0.05,
verbose = TRUE,
OrgDb = organism,
pAdjustMethod = "none",
eps = 0.0)
dotplot(gse, showCategory=10, split=".sign") + facet_grid(.~.sign)